Multielectron corrections in molecular high-order harmonic generation for different 

formulations of the strong-field approximation 
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We make a detailed assessment of which form of the dipole operator to use in calculating high order 
harmonic generation within the framework of the strong field approximation, and look specifically 
at the role the form plays in the inclusion of multielectron effects perturbatively with regard to the 
contributions of the highest occupied molecular orbital. We focus on how these corrections affect the 
high-order harmonic spectra from aligned homonuclear and heteronuclear molecules, exemplified by 
N2 and CO, respectively, which are isoelectronic. We find that the velocity form incorrectly finds 
zero static dipole moment in heteronuclear molecules. In contrast, the length form of the dipole 
operator leads to the physically expected non-vanishing expectation value for the dipole operator 
in this case. Furthermore, the so called "overlap" integrals, in which the dipole matrix element is 
computed using wavefunctions at different centers in the molecule, are prominent in the first-order 
multielectron corrections for the velocity form, and should not be ignored. Finally, inclusion of 
the multielectron corrections has very little effect on the spectrum. This suggests that relaxation, 
excitation and the dynamic motion of the core are important in order to describe multielectron 
effects in molecular high-order high harmonic generation. 



I. INTRODUCTION 

High harmonic generation can be used for attosecond 
imaging of molecules Q- This is a consequence of the 
fact that the structure of an aligned molecule can be in- 
ferred from quantum-interference patterns, which appear 
in the high-harmonic spectra. Understood in terms of 
the three-step model, in which an electron reaches the 
continuum by tunneling or multiphoton ionization, is ac- 
celerated by the field and subsequently recombines with 
a bound state of its parent molecule, emitting a high- 
order harmonic photon |l|, this phenomenon occurs due 
to the electron returning to different spatially separated 
centers. The simplest scenario is high-harmonic genera- 
tion in a diatomic molecule, for which, in principle, re- 
combination and thus high-harmonic emission can take 
place from either atomic center. This is a microscopic 
equivalent to the double slit experiment Q. 

How to accurately describe the whole process of high 
harmonic generation (HHG) in molecules is still an open 
problem. There are two main theoretical approaches that 
can be used. Firstly one can attempt a fully numerical so- 
lution by using density functional theory, which gives ac- 
curate results but inhibits a more complete understand- 
ing of the process [J] . Secondly, one can employ the tra- 
ditional atomic three step model, within the strong-field 
approximation [l|, and calculate the molecular bound- 
state wavefunctions using a quantum chemistry code [6j . 
In this framework, the bound states of the molecule are 
mainly incorporated in the ionization and recombination 
steps, in the form of prefactors. The latter approach al- 
lows a better physical understanding of the problem. The 
price one pays, however, is that, in order to be able to 
treat HHG within a semi-analytical framework, several 
approximations must be made. 

It is well known that the harmonic spectrum is propor- 



tional to the modulus squared of the Fourier transform 
of the dipole acceleration. By using the Ehrenfest theo- 
rem the dipole acceleration can be replaced by the dipole 
velocity or the dipole length, with no loss of accuracy so 
long as the integration is over all time and no physical ap- 
proximations are made upon all wavefunctions involved. 
However, when the SFA is used, this is no longer the case, 
and replacement of the dipole operator gives rise to incon- 
sistencies in the SFA dipole matrix elements, which vary 
depending on the form of the operator used. The most 
well known of these is the lack of translational invariance 
when using the length form of the dipole operator. In the 
single-active electron approximation, it has been shown 
that the velocity form gives the most reliable results [HJ 
and is easier to use than the acceleration form. 

In calculating the SFA prefactors the simplest ap- 
proach is to use the highest occupied molecular orbital 
(HOMO) j5|. However, the problems that arise in this 
case are 1) which orbitals have the predominant influence 
on high harmonic generation and 2) how to accurately in- 
clude the multielectron effects of electron exchange and 
correlation. 

The first problem has to some extent been addressed 
in Ref. [2[ by noting that the ionization potential of 
different orbitals varies with the rotation of the molecule 
and hence the corresponding cutoffs in the high harmonic 
spectra will be different. Furthermore, the symmetry of 
the orbital will strongly affect the contribution of specific 
orbitals to high-order harmonic spectra. This is partic- 
ularly true as the alignment angle of the molecule with 
respect to the laser-field polarization varies. In previ- 
ous work [8| we also investigated high-order harmonic 
generation beyond the single-active orbital approxima- 
tion. In particular, we addressed the quantum interfer- 
ence between the HOMO and the HOMO-1 of N 2 and 
the HOMO and the lowest unoccupied molecular orbital 



(LUMO) of N2 • For the latter case, we have also ob- 
served that electron recombination to different orbitals 
can be mapped into different cutoffs in the HHG spectra. 
Our investigations, however, did not incorporate electron 
correlation. In fact, we have either dealt with simplified 
multielectron models, in which the time evolution of the 
optically active electrons has been decoupled, or with a 
coherent superposition of one-electron states. 

The second problem is what we propose to investigate 
in this paper and in doing so to consider the influence 
that these multielectron effects have on high harmonic 
generation. In particular, we wish to incorporate such ef- 
fects perturbatively around the contributions of the high- 
est occupied molecular orbital. Although correlation in 
molecular high-order harmonic generation has been ad- 
dressed in Refs. @, la], utilizing a multielectron version 
of the strong-field approximation, such studies have been 
performed numerically to a great extent. Including the 
correlation analytically in the dipole matrix elements fa- 
cilitates more insight into the effects that correlation has 
on the high harmonic spectra. In Ref. [7|, many-body 
perturbation theory was used to include electron corre- 
lation as corrections to the standard three step model. 
Therein, however, only the response of single atoms has 
been addressed. 

From the viewpoint of which form of the dipole oper- 
ator to use, we are particularly interested in how such 
corrections are affected by the different formulations, es- 
pecially in molecular targets. Indeed, even in the single- 
active electron approximation, the discrepant results ob- 
tained for each dipole form have raised considerable de- 
bate. By extending the debate to multielectron effects 
in molecules, more insight into which operator is most 
accurate for harmonic generation, and indeed all strong 
field phenomena, is gained. The use of molecules is ben- 
eficial because the two-center interference pattern gives 
another indicator of how the formulation is performing. 

In order to make such an assessment, we will em- 
ploy many-body perturbation theory along the lines of 
Ref. pi and the dipole moment in its length and velocity 
forms. For that purpose, we will consider two isoelec- 
tronic molecules: N 2 and CO. The HOMO, LUMO and 
HOMO-1 in both molecules have similar geometry. The 
main difference between them is that N2 is a homonu- 
clear molecule, while CO is heteronuclear. This means 
that, for the latter molecule, there is an intrinsic static 
dipole moment. According to Ref. p}, one expects the 
multielectron corrections to be more prominent in this 
case. 

This work is organized as follows. In Sec.|ITl we provide 
the theoretical framework which we will employ in order 
to compute the high-order harmonic spectra. We then 
calculate the transition amplitudes, in different formula- 
tions (Sec. Ill A"|) . within the strong field approximation, 
and calculate the multielectron corrections (Sec. IIIBD . 
with particular focus on the form of the dipole operator. 
Subsequently, in Sec. Mil these expressions are employed 
to compute the high-harmonic spectra. Finally, the pa- 



per is summarized in Sec. IIVI Atomic units are used 
throughout. 



II. THEORY 

A. Transition amplitude 

The SFA transition amplitude in the particular formu- 
lation of Ref. [l] reads 

/'OO ft /> 

dt / dt' / d 3 fc<oc(k + A(t))a ion (k + A(t')) 
-00/ — OO J 

exp[iS(t,t',Cl,k)] + c.c, (1) 

with the action 

1 /•* 

S(t,t',n,k) = -~ / [k + A{T)] 2 dT-e Q (t-t') + nt (2) 

2 Jt' 

and the prefactors a rec (k + A(i)) = (k + A(i)| d • e x |\f) 
and a ion (k + A(f )) = (k + A(t')| E(f ) • r |*) . In Eqs. 
((T|) and ©, d, e x , eo, and f2 denote the dipole operator, 
the laser-polarization vector, the ionization potential of 
the molecule in question, and the harmonic frequency, re- 
spectively, and I \P) gives the bound state with which the 
electron recombines, or from which it tunnels. Eq. (TTJ) is 
given in the length-gauge formulation of the SFA. We will 
employ this gauge throughout. In the standard velocity- 
gauge formulation for the molecular SFA, the two-center 
interference patterns vanish [ll|, Il5j . To first approx- 
imation, we assume that there is one active electron, 
which tunnels from the highest occupied molecular or- 
bital, reaches the continuum and recombines. Many- 
body corrections due to the presence of the other elec- 
trons will only be incorporated in the prefactors. 

The above-stated expression describes a physical pro- 
cess in which an electron tunnels from its parent molecule 
at an instant t' , and propagates in the continuum from 
a time t' to a subsequent time t, with the intermedi- 
ate momentum k. At t, the electron recombines with 
its parent ion, generating a high-energy photon of fre- 
quency f2. These steps are implicit in the action ([2]). They 
can also be directly identified in the saddle-point equa- 
tions, obtained from the values of t, t' and k which render 
the action stationary. This implies dS(t,t',£l,\iL)/dt' = 
dS(t, t', Q, k)/dt = and dS(t, t', ft, k)/fik = 0. 

These saddle-point equations read 



[k + A(t')Y 



+ e = 0, 



/ dr [k + A(t)] = 0. 
Jt' 



and 



[k + A(i)] 



e = fi. 



(3) 



(4) 



(5) 



Physically, Eq. ([3]) gives the kinetic energy of the electron 
during tunneling. One should note that this equation has 
no real solution, so that the tunneling time t' is complex. 
This is a consequence of the fact that tunneling has no 
classical counterpart. Eq. Q constrains the momentum 
of the electron, so that it may return to the geometrical 
center of its parent molecule. Finally, Eq. ([5]) gives the 
conservation of energy upon recombination, in which the 
kinetic energy of the returning electron is converted in 
a high-frequency photon. The stationary-phase method 
will be employed to compute the spectra in this work. 
For more details we refer to [12]. 



B. Many-body corrections 

We will now briefly discuss the many-body corrections 
derived in Ref. T'j . They have been obtained by starting 
from a many-body dipole operator in the recombination 
and ionization dipole matrix elements a rcc (k + A(t)) and 
aj on (k + A(i')), and applying Moller Plesset [13J pertur- 
bation theory. 

a„(k) = af (k) + a«(k) + a^(k) + 0(H 2 ), (6) 

where the dipole matrix element associated with the 
HOMO is given by 

af(k) = <k|d|Vo). (7) 

The first-order corrections read 

a«(k) = £)(d«<k|^) - d i0 (k|^}) (8) 



and 



af(k) = ^)EE^ 
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where 



d,^ = (^|d|-0 M ), 



(9) 



(10) 



v aO[ij] = (lpi,1pj\v\lpa,1po) ~ (^j,1pi\v\lp a ,ipo) , (11) 



and 



v ab[0j] = {^0,^j\v\lj)a,i>b) ~ (lpj,1po\v\'>pa,'>Pb) ■ (12) 

In the above-stated equations, v denotes the electron- 
electron correlation, e a the orbital energies, and the index 
•q, with 77 = (rec, ion), refers to recombination and ioniza- 
tion. The indices a, b represent unoccupied orbitals and 
the indices i,j represent occupied orbitals. The index 
is related to the orbital around which the corrections are 



inserted. In the specific scenario addressed in this paper, 
this is the highest occupied molecular orbital (HOMO). 
The indices v, [i in the dipole matrix element (fTU)) are 
general. In the first-order corrections ©, they may 
relate to the expectation value of the dipole operator 
{y = jj, = i), or to the dipole coupling between the 
HOMO and an occupied bound state [y = i and (x = 0). 
In the framework of the second-order corrections ©, 
they relate to the couplings occupied and unoccupied 
bound states (y = a and /i = j or v — i and fi = b). Phys- 
ically, v a o[ij] an d v ab[0j] gi ye the exchange terms, and are 
dependent on the form of the electron-electron interac- 
tion. In the corrections (|8]) and ([9]), the terms (k|?/v) will 
determine the two-center interference conditions, while 
the remaining terms will mainly act as weights. 

These corrections will also behave differently with re- 
gard to homo nuclear and heteronuclear molecules. For 
the former type of molecules, the orbitals exhibit definite 
parities. Hence, the dipole matrix elements d^ will only 
couple gerade and ungerade states. As a direct conse- 
quence, the first term in ((5J will vanish. Physically, this 
is expected, as there is no static dipole moment for a 
homonuclear molecule. 

In contrast, for heteronuclear molecules, the matrix 
elements d y/J will in principle couple all states involved. 
The first term in Eq. (|SJ) will be non- vanishing, and, 
physically, is the contribution of the static dipole moment 
of the molecule in question. This term, however, will only 
lead to quantitative changes in the dipole matrix element. 
In fact, an inspection of Eqs. ([7]) and (|8j) shows that it 
will yield the same two-center interference condition as 
the one- particle prefactor a v (k) . 

In the above-mentioned framework, all the structure of 
the molecule is embedded in the recombination prefactor 
a rcc (k + A(t)). This prefactor can be written in different 
forms, which will lead to different results. Explicitly, 
the length, velocity and acceleration forms of the dipole 
operator read 



and 



d (0 



d^ = -k 



d (a) = -vm 



(13) 



(14) 



(15) 



respectively, where the hats denote operators. The pref- 
actor aj on (k + A(£')) will mainly influence the overall 
harmonic intensity, as it is associated with ionization. 
In this work, we will consider the length and the velocity 
forms of the dipole operator. In the single-active-electron 
framework, the former and the latter lead to the worst 
and best description of the two-center interference pat- 
terns [101 ] . 



C. Molecular orbitals and interference condition 

In calculating the molecular orbitals we use a Linear 
Combination of Atomic Orbitals (LCAO) approximation 
along with the Born Oppenheimer approximation. In 
this framework, the molecular orbital wavefunction is 



Mr) = E c S $ S( r + R / 2 ) 



+(-1) 



l a -r 



■;0 



c aj $^(r-R/2) (16) 



where R is the internuclear separation, l a and m a the 
orbital and magnetic quantum number respectively. The 
coefficients c a -, with £ = L, R, form the linear com- 
bination of atomic orbitals which are extracted from 
GAMESS-UK [6j, and the indices (L) and (R) refer to 
the left or to the right ion, respectively. The internu- 
clear axis is taken to be in the z direction, and the 
laser-field polarization is chosen along the radial coor- 
dinate. For homonuclear molecules, c a = c^ ■ = c QJ - 

and & a • = $ a l = ®a,j- The parameter A Q can be varied 
depending on the symmetry. For homonuclear molecules, 
A a =|r7i Q | and Aq,=|771q,|-(-1 for gerade and ungerade sym- 
metry, respectively. 

The wavefunctions themselves are then expanded as 
Gaussian type orbitals with a 6-31 basis set, with 



$®.(r) 

a, 7 V / 



EC- ■'• 



(r, 



,-C 



(17) 



Hence, for a linearly polarized field 9k is equal to the 
alignment angle 9l between the molecule and the field. 
The above-stated equations have been derived under the 
assumption that only s and p states will be employed in 
order to build the wavefunctions employed in this work. 
For more general expressions see Rcf . [8] . 

In Fig. 1, we exhibit the HOMO for and CO and N 2 
obtained with GAMESS-UK [Figs. [U (a) and (b), respec- 
tively]. The position-space wavefunctions exhibit a cen- 
tral maxima and two minima located at the ion for both 
species. The main difference is that, for CO, there is a 
strong bias towards the Carbon atom, while for N 2 the 
wavefunction is symmetric. This is in agreement with the 
standard molecular-physics literature. In contrast, we 
have found that both momentum-space wavefunctions, 
given in Figs. |TJ(c) and (d), are symmetric with respect 
to (pxtPz) ~~> (~Px,—Pz)- In comparison to N 2 , CO ex- 
hibits a much deeper maximum at vanishing momenta. 
Such maxima and minima may be determined by writ- 
ing the exponents in Eq. (fiT))) in terms of trigonometric 
functions. In this case, one obtains 



^■(k)=£4 a, cos 



k R 



iC sin 



k R 



. (21) 



with 

C { ± ] = (-1)'-— +A "cS<O k ) ±cg*S(k). ( 22 ) 
Calling ■& = arctan[iCi a) /Ci a) ],we find 



For the a, ir x and ir y orbitals, r x = z, x and y, respec- 
tively. The index £ is related to the left or right ion, and 
b v ,j and Q ■ form the contraction coefficients and expo- 
nents, also obtained from GAMESS-UK. Note that these 
coefficients are real. 

The momentum-space wavefunctions, which will be 
used in this work to compute the prefactors, are given 
by 



^•(k) 



^exp[*k.^] C g$g(k) 






(18) 



R, 



''' '" v "expHk~]c^$|3(k), 



where, similarly to the position-space situation 

*S(k)=E 6 S^S( k ). ( 1Q ) 



with 



r 3/2 



(ih 



(<8) 



I7 ^ex P [-fc7(4C^)]. 



( 2 °) 

Therein, /3 = z, = X and (3 — y for the a, tt x and it y or- 
bitals, respectively. The return condition Q guarantees 
that the momentum k and the external field are collinear. 



^•(k) 



E 



>(«) 



& ] ) sin^ + k-R/2]. (23) 



rm. 



Eq. (J23J) exhibit minima for -d + k • R/2 

Note, however, that the coefficients C±' defined in 
Eq. (|22|) depend on the wavefunctions at the left and 
right ions. Since these wavefunctions themselves de- 
pend on the momentum k, one expects the two-center 
patterns to be blurred for heteronuclear molecules. In 



contrast, for homonuclear molecules, c 
$l L ](k) = $g(k). 



(£) 



c^j and 

>„y(k) = $q, (k). This implies that the momentum 
dependence in the argument •& cancels out and that the 
interference condition in Refs. [a, [|| is recovered. In 
this case, sharp interference fringes are expected to be 
present. 

The above-stated condition does not only lead to inter- 
ference fringes in the bound-state momentum wavefunc- 
tions but also in the high-harmonic spectra. This is due 
to the fact that the dipole matrix elements depend on the 
wavefunctions (fT9|) . In this latter case, k — > k + A(£) in 
Eqs. flm-flU). 



D. Dipole matrix elements 

In this section, we will provide explicit expressions for 
the dipole matrix elements, depending on whether the 




with 



;(k) 



R, 



^exp[zk.|]cS$g(k) 



(26) 



-(-!)'" 



J-R) 



R, 



c^exp[-ik~]*^(k) 



'"•■The first two terms in Eq. (|2"5|) are mainly oscillating, 

! ™and give the two-center interference maxima and minima. 

"°The third term in Eq. (|25|) increases with the internu- 

"""clear distance, and has the main effect of blurring such 

patterns. In fact, the existence of such term has caused 

considerable debate and has been attributed to the lack 

of orthogonality between the continuum and the bound- 

'""state wavefunctions in the SFA [14J . Similar terms are 

""also present in the first-order corrections, to be discussed 

slater. 

In the velocity form, the single-electron dipole matrix 
element reads 



,(") 



R 



X)C>exp[ik.=]$W(k) 



(27) 



FIG. 1: Highest occupied molecular orbitals for CO [panels 
(a) and (c)] and N2 [panels (b) and (b)], in the position and 
momentum space (upper and lower panels, respectively). 



dipole operator is given in its length or velocity form. 
We will restrict ourselves to the first-order corrections to 
the single-active-electron prefactor, as they are expected 
to be dominant. 



+(-l) I -- m - +A '"cgkexpHk.|]#2(k). 

In contrast to the length-form scenario, there is no term 
linearly growing with the internuclear separation. In- 
deed, previous studies for simplified models for homonu- 
clear molecules indicate that this is the preferred form of 
the dipole operator, in order to obtain the correct inter- 
ference patterns [lOt ] . 



2. Corrections a^ (k) 



1. Single-active electron approximation 

We will commence with the single-active-electron 
dipole prefactor ([7]), which is employed in the standard 
strong-field approximation. In the length form, this pref- 
actor reads 



a 1 (0) (k)=- l a k ^(k), 



(24) 



where ?pj (k) is given by Eq. (|19l) . If explicitly written in 
terms of the wavefunctions centered at the left and right 
ions, 



a( 0) (k) 



R, 



X C gexp[ l k.-]a k <j(k) 



+ (-1) 

R 

1 



l a -m a +\ a (R) 



+ -&(*■), 



exp[-ik-|]^g(k) 

(25) 



We will now write down the explicit forms of Eq. (j8]) 
for molecular orbitals in the LCAO approximation. Ac- 
cording to Ref. [7] , this term is responsible for the domi- 
nant corrections to the single-electron prefactor, and de- 
pends on the momentum space wavefunction (|19p . and 
on the matrix element (|1U[) . The former is independent 
on the form of the dipole operator, while the latter is 
explicitly given by the sum 



<W = Y^ Y Yl E £'>dla, h)Z£ft: 

<x,i<*P,h ?',« 



(28) 



where 



r a,0 



= d 



-(€') 






lK ( P(r + (-irn/2) 



xd(r)$f(v+(-l) b n/2). 



(29) 



In Eq. ([51) and (gSJ, the indices £',£ are related to 
the right or left ions, the exponents a, b are 1 and 
for £',£ = R and £',£ = L, respectively, and d(r) de- 
notes the dipole operator in the position representa- 
tion. In its length and velocity forms, d^(r) = r and 



d(")(r) = iV r .The coefScients E^i^(l a ,lp) are such that 

S£,i(J o ,J/0 = l,S£,K(«a,^) = (-i)'" +A/ »,Sji jL (i a ,/ / ,) = 
(_!)I-+A n and Z R<R (l a ,lp) = (_l)I/.+Afl+I-+Aa < 

In Eq. (|29[) . one may recognize two types of contribu- 
tions. The direct integrals 



x €,€ 



dfr 



,(£) 






(0 



$Hf?(r±R/2)d(r)$^(r±R/2), 

(30) 

which contains one-electron wavefunctions centered at 
the same ion, and the so-called "overlap" integrals 



T a,p 



<Tr 



M') 






$^J(r±R/2)d(r)$^(r T R/2), 

(31) 

for £' 7^ £, with one-electron wavefunctions localized at 
different ions. For the moment ignoring the overlap inte- 
grals, 









(32) 



I /^y^+A^+i Q +A 



,(fl) 



c P, f j,- L R,R- 



The explicit expressions for the overlap integrals are pro- 
vided in Appendix 1. 

We will commence by providing expressions for (|28[) us- 
ing the velocity form. For that purpose, one must solve 
the integrals (|2T))) . Depending on whether the dipole op- 
erator couples only a or n orbitals, or a a orbital with 
a 7r orbital, these integrals will be different. If only a 
orbitals are involved, 



2gW)=£- 



7T&W 



(33) 



x 



[ifiHlfi +i a -i)~ xf^Hh + la + i)} , 
with £ = i? or L, and 

•f(o = J [i +(-!)'] (^S+®" va "' /arn + / 



(34) 
The same expression is encountered if only ir x or ir y are 

coupled, i.e., l£'f(a,a) = l^f(n x ,n x ) = 2? ^ (%,%)■ 
The above-stated expression is only non- vanishing if Z^ + 
Z Q + 1 is even. In our framework, this implies that only s 
and p states are mixed. If, on the other hand, a a orbital 
is coupled to a n x or a 7Tj, orbital, this integral reads 






(o'.TTx) 






1 (c (i) + 6P ) 



1/2 



x(6(/ Q ,^)-2CJ^ Q + l,Z /3 )),(35) 



with x = x, y. Thereby, 

A(i a ,lp) =FQa)Hb), 



(36) 



B(l a ,lp)=T{l a ){lpF(lp-l) 



K^Hh 



i 



>}• 



(37) 



and F{lp) is defined according to Eq. (|3"4"|) . Eq. (J3~5j) 
also holds for the dipole matric elements coupling the 
tt x and TT y orbitals, i..e, I?'f (<t,7T x ) = Iff (n x ,TT y ) = 

T?f(ir y , tt x ). From Eqs. (f36f and (|37| we see that ,4(Z a + 
1, //j) ^ if Z Q = and ^ = 1, and that B(l a , Ip) ^ if 
Z Q = 1 and Ip = 0. Hence, once more, only s and p states 
are coupled. 

All these contributions are then included in the cor- 
rections ©. In the specific case of a homonuclear 
molecule, an inspection of Eqs. (|32 |) - (|37p shows that 
only orbitals of different parities are coupled. In fact, 

it c a ^v = c a ,i/ and Co 



Ca L Eq. (|3"2")l is only nonvan- 



ishing if Ip + \p + l a + A a is even. Since, as discussed 
above, If f are only nonvanishing if l a and Ip have values 
corresponding to different parities, this must also hold 
for \p and A Q . In contrast, for heteronuclear molecules, 
these coefficients are different and in principle all terms 
contribute. 

We will now state the length-form dipole matrix ele- 
ment. We will mainly consider the integrals Iff . If only 
a orbitals are coupled, these integrals read 



Zlf(cr,cr) 






TTfcW 
3, v 3 >A« 

3,3' v»3>* s 3 ^' 

J 7 (l a + 1/3 



(38) 



R 



1)T ^(L + lfi) 



The same expression holds for transitions involving only 
the n x or the ir y orbitals, i.e., Iff (a, a) — Iff(n x , ir x ) = 

If the dipole operator couples a with tt orbitals, such 
integrals are given by 
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with x = a;, V- Eq. (J40I) also holds for couplings between 
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different 7r orbitals, i.e., I?'; (a, 7T S 

I iA (7Tx,7r y )- 

In the above-stated expressions, one may identify two 
distinct behaviors. The terms ^(la + lp + l), A(l a + l,lp) 
and A(l a ,lp + 1) are non- vanishing only if l a and lp 
have values corresponding to different parities, i.e., in our 
framework, it couples s and p states. The terms IF(l a +lp) 
and A(l a ,lp), on the other hand, couple states of the 
same parity (either s or p). For homonuclear molecules, 
these terms will vanish, as the dipole operator will only 
couple states of different parities, i.e., Ip + l a + 1 must be 
even in Eq. (|3"2"j) . For heteronuclear molecules, in con- 
trast, they may in principle be present as the coefficients 
at each center of the molecule are different. 



E. Inconsistencies in the dipole matrix elements 

We will now discuss some of the subtleties that arise in 
the evaluation of our multielectron corrections. Firstly, 
the overlap integrals, in Eq. pip are often presumed to 
be very small in comparison to their direct counterpart, 
but this is not always the case. This can be seen by using 
the Gaussian theorem, 
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Here £ Q and Qp correspond to our exponents in Eq. (|17|l . 

At first sight it would appear that these contributions 
are much smaller than those from the direct integrals 
(150)) . This does not hold, however, for exponents that 
are less than unity, of which there are many for various 
different atomic orbitals. It is therefore not astute to 
ignore the overlap integrals in this instance. 

In fact, we have verified that in most cases the over- 
lap integrals (l3"Tj) are comparable to the integrals (f3"tJ)) at 
the same center. More importantly, we have also found 
that the velocity form gives zero static dipole moment 
for heteronuclear molecules when using a Gaussian basis 
set. This is an unphysical result, which has been verified 
for the direct integrals by using the properties of Gamma 
functions to analyze the term in brackets in Eq. (|33|) . Ex- 
plicitly, it can be shown that this term is proportional to 
[Qlj3 — (jl a /(('j + Cj)]- Since only s and p states can be 
coupled, each pair (l a ,lp) = (0,1) and (l a ,lp) = (1,0) 
will lead to symmetric contributions which will cancel 
out. We have also verified numerically that the overlap 
integrals are vanishingly small in this case. 

Apart from that, inclusion of the overlap integrals gives 
rise to terms that depend on R, and thus removes the 
translational invariance from the velocity form. This 
problem is also encountered in the single-active electron 
matrix element when using the length form. In contrast 
to the one-electron length-form scenario, however, these 
terms will vanish at large internuclear distances due to 
the Gaussian exponentials (1421) . (see a discussion in [14| 
for their counterpart in the single-active electron frame- 
work) 

In contrast, the length form predicts non-zero static 
dipole moment in heteronuclear molecules, which makes 
sense physically. However, for the heteronuclear case the 
R-dependent terms that arise in the direct integrals cou- 
ple states of the same parity as seen in Eq. (IMj) and 



Eq. (1401) . This should not occur for an odd operator 
such as the dipole, which is an unphysical result. This 
is another reason to suggest that these terms should be 
neglected. In the results that follows, this term will be 
neglected unless otherwise stated. 



III. HIGH-ORDER HARMONIC SPECTRA 

Below we compare the harmonic spectra predicted by 
different forms of the dipole operator. We will place par- 
ticular emphasis on the comparison between homonuclear 
and heteronuclear molecules. The former and the lat- 
ter will be exemplified by N2 and CO, respectively. We 
will consider the full, three-dimensional dynamics of the 
problem. For that purpose, one must integrate over the 
azimuthal coordinate, <j) p , such that, 
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where M is the overall transition amplitude. This al- 
lows one to consider the degeneracy of the 7r orbitals [8] . 
The corrections are incorporated in the ionization and 
recombination prefactors. The former gives the partic- 
ular weight for tunneling to take place, and the latter 
gives the shape of the high-order harmonic spectra and 
the two-center fringes. The explicit expressions for the 
dipole matrix elements are provided in Appendix 2. 

We will now start by discussing the high-harmonic 
spectra of diatomic nitrogen (N2) computed including 
the single-atom prefactor and the first-order corrections. 
In Fig. [21 we display such spectra, in the velocity and 
length form (panels (a) and (b), respectively), together 
with the corrections alone (panels (c) and (d), respec- 
tively). In the length- form prefactor a a (k), we omitted 
the term growing linearly with the internuclear distance 
R. In the corrections, we included both the direct and 
the overlap integrals. 

In the figure, we notice a clear two-center mini- 
mum, which is due to the interference of high-harmonic 
emission at the two spatially separated centers in the 
molecule. This minimum varies with the alignment angle 
6l of the molecule with respect to the laser-field polariza- 
tion. In the upper panels, its energy position is charac- 
teristic of the 3a g orbital |8|, |T5[ . A comparison of panels 
(a) and (b) shows that the two-center minimum is slightly 
displaced, depending on the form of the dipole operator. 
For instance, for vanishing alignment angle, in the ve- 
locity form this minimum is near O = 21a; , while in the 
length form it is at a slightly higher energy ($7 = 25w). 
This is due to the fact that the s — p mixing, which deter- 
mines the energy position of this minimum, is different 
for each case. The velocity form slightly favors the contri- 
butions from the p states, compared to the length form. 
For a discussion of the s — p mixing in a slightly different 
context, namely of different basis sets, and its influence 
on the high-order harmonic spectra, see |8j. Here, how- 
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M0 FIG. 3: High-order harmonic spectra for N2 in a driving field 
jvith the same parameters as in the previous figure, and an 
alignment angle 9l = 7r/4. The figure shows the overall spec- 
'tra together with the first-order corrections taking into ac- 
count only the direct or overlap integrals. We consider both 

. 270 £he ionization and recombination prefactors, or the recom- 
bination prefactor only (orange and red lines in the figure). 
Panels (a) and (b) give the velocity and length forms of the 
dipole operator, respectively. 



FIG. 2: High-order harmonic spectra for N2 subject to a lin- 
early polarized laser field of frequency uj — 0.057 a.u. and 
intensity / = 4 x 10 14 W/cm , as a function of the alignment 
angle 6l between the molecule and the field. In the figure we 
consider the corrections in both the ionization and the recom- 
bination prefactors. Panels (a) and (c) have been computed 
employing the length form, while panels (b) and (d) have 
been calculated using the velocity form. The upper panels 
display the 3<r 9 spectrum in which the first-order corrections 
have been incorporated. In the lower panels, we depict the 
yield obtained with the first order corrections in the length 
and the velocity form. The bound-state energy of the HOMO 
and the equilibrium internuclear distance have been taken 



from GAMESS-UK. Explicitly, e 



(N 2 ) 



0.63485797 a.u and 



RS " = 2.068 a.u. The contour plots have been normalized 
with regard to the maximum yield in each panel. 



ever, this bias is introduced by the different forms of the 
dipole operator. 

The remaining panels exhibit the spectra obtained em- 
ploying the first-order corrections only. Such contribu- 
tions exhibit a well-defined minimum, whose position is 
at a quite different energy from that observed in the over- 
all prefactor. We have verified that this pattern can be 
related to the contributions of the a u orbitals. The con- 
tributions of the 7T„ orbitals are strongly suppressed in 
this parameter region due to to the presence of a nodal 
plane along the internuclear axis. 

Quantitatively, the yield in panels (c) and (d) is sev- 
eral orders of magnitude smaller than those in the upper 
panels. This shows that, in the present framework, it 
suffices to consider the zeroth order contributions. This 
behavior occurs both for the length and the velocity form 
of the dipole operator, and is due to the fact that the cor- 
rections are considerably smaller both in the ionization 



and recombination prefactor. The above-mentioned sup- 
pression of ionization for the it orbitals also contributes 
to such results. 

The influence of the ionization prefactor can be seen in 
Fig. [3] Therein, we present spectra computed for a fixed 
alignment angle assuming aj on (k + A (to)) to be constant 
and equal to unity. Physically, this implies that any in- 
fluence of the orbital symmetry on tunneling ionization, 
such as the presence or absence of nodal planes, has been 
removed, both for the zeroth-order prefactor and the cor- 
rections. For comparison, we are also providing the full 
spectra, with both the ionization and recombination pref- 
actors. 

The figure shows that, while the overall spectra 
changes relatively little, or may even increase if this pref- 
actor is removed (see orange and black lines in Fig.[3j(a)), 
the corrections increase in several orders of magnitude. 
This can be attributed to the following. In the parameter 
range of interest, the zeroth order ionization prefactor, 
which is the main influence in the overall spectra, is of 
the order of unity or even slightly larger for a 2>a g orbital. 
This is due to the absence of nodal planes for either par- 
allel or perpendicular alignment in this particular case. 
In contrast, in the corrections, mainly a u or tt u orbitals 
are involved. The former orbitals exhibit nodal planes 
for alignment angles 9l = 90° and the latter for 9l = 0. 
This will cause an overall suppression in tunnel ioniza- 
tion for a wide range of angles, when such orbitals are 
summed over. Furthermore, the dipole matrix elements 
d^ in the spectra are quite small and will also contribute 
to the above-mentioned suppression. 

Another noteworthy feature also presented in Fig. [3] is 
the influence of the so-called overlap integrals, in which 
the dipole operator couples wavefunctions at different 
centers in the molecule, in the multielectron corrections. 
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FIG. 4: High-order harmonic spectra for CO in a driving 
field with the same parameters as in the previous figure, as 
a function of the alignment angle. Panels (a) and (c) exhibit 
the full spectra and the corrections in the length form, while 
panels (b) and (d) give their counterparts in the velocity form. 
In the length form, the linearly growing terms with regard to 
the internuclear distance have been neglected. 



In general, these integrals are assumed to be very small in 
comparison with the direct integrals, in which wave func- 
tions at the same ion are taken. Hence, in many cases 
the overlap integrals are neglected. We have verified, 
however, that this depends strongly on the form of the 
dipole operator. For the velocity-form dipole, the spec- 
tra obtained from the overlap integrals are comparable to 
those obtained from the direct integrals only [Fig. [31(a)]. 
Hence, both contributions must be included. In contrast, 
if the length form is considered, the contributions from 
the overlap integrals are considerably smaller than their 
direct- integral counterparts [Fig. |3J(b)]. Therefore, the 
former may be neglected. 

Next, we will discuss the spectra encountered for CO. 
The overall spectra are depicted in Fig. HI together with 
the first-order corrections (upper and lower panels, re- 
spectively). As an overall feature, the spectra exhibit no 
perceivable two-center interference minimum, as CO is a 
heteronuclear molecule. This means that, upon recom- 
bination, there will be more high-order harmonic emis- 
sion from a specific spatial region, where the orbital with 
which the electron recombines, i.e., the HOMO, is more 
localized. This occurs both for the velocity and length 
forms (Figs. |U (a) an d (b), respectively). 

One notices, however, that the decrease in the over- 
all yield as the molecular alignment angle approaches 
9l = 90° remains, due to the fact that the HOMO is 
a a orbital, and thus highly localized along the internu- 
clear axis. As the laser-field polarization direction moves 



away from this axis, recombination, and hence high-order 
harmonic emission becomes less and less probable. Due, 
however, to the heteronuclear nature of this molecule, 
this suppression occurs at a different angle compared to 
N 2 . 

Once more, we find that the contributions of the first- 
order multielectron corrections to the spectra are several 
orders of magnitude smaller than the zeroth-order terms. 
These contributions are depicted in Fig.|4j(c) and (d), for 
the length and velocity form, respectively. For such cor- 
rections, the interference minimum due to high-harmonic 
emission at different centers has been washed out. This 
minimum was very clear for N2 (see Figs. [2] (c) and (d) 
for comparison). There is, however, a very pronounced 
minimum near 8^ = 50° for the velocity-form correc- 
tions, which is independent of the harmonic frequency 
(see Fig. Hl(d)). This minimum is due to the geometry 
of the three lower-lying a orbitals, which influence the 
corrections. This suggests that the imprints on the high- 
order harmonic spectra due to purely geometric features, 
such as nodal planes or minima in the orbitals, are more 
robust than those related to quantum interference be- 
tween harmonics emitted as spatially separated centers. 
The length-form corrections, in contrast, do not exhibit 
any such feature. In fact, we have verified that they are 
dominated by the static dipole moment, which behaves 
in the same way as the HOMO. 

Below, in Fig. [5] we perform a more detailed analysis of 
the above-mentioned issue. For that purpose, we consider 
the contributions of individual orbitals to the corrections, 
and their influence on the harmonic spectra, for fixed 
harmonic order. We consider a mid-plateau harmonic, 
namely fl = 45w. These results are displayed for both N2 
and CO, as functions of the alignment angle (upper and 
lower panels in the figure, respectively). As an overall 
feature, we observe that, for a homonuclear molecule, 
all contributions are symmetric upon 9 — > 8 + tt (see 
panels(a) and (b)). Physically, this is expected, as such 
molecules exhibit inversion symmetry and a rotation of 
tt should not alter the physics of the problem. Clearly, 
this no longer holds for a heteronuclear molecule. In this 
latter case, however, the contribution of S n {6i) to the 
yield is the mirror image of S u (9l + 7r). This is clear, as 
in this case the right and the left ions are reversed. 

The individual contributions provide useful informa- 
tion about the geometry of the orbitals involved in the 
corrections. This can be seen in all cases, regardless of 
whether the molecule is heteronuclear or homonuclear. 
For instance, the contributions from the 7r orbitals, shown 
as the orange symbols in the figure, decay very abruptly 
for 8l = nir. This is due to the fact that such orbitals ex- 
hibit nodal planes in the case of N2 and a feature similar 
to nodal planes in the case of CO. The main difference 
between the homonuclear and the heteronuclear case is 
the two-center interference minimum, which is absent in 
the heteronuclear case. Furthermore, the shapes of such 
distributions do not change with the form of the dipole 
operator. Any changes are mostly quantitative (for a 
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FIG. 5: Contributions of individual orbitals to the first-order 
corrections in the 45 t h harmonic of a field with the same pa- 
rameters as in the previous figures, as functions of the align- 
ment angle 9l of the molecule relative to the laser-field po- 
larization. The thick black lines in the panel give the overall 
contributions of all a orbitals and the orange stars the con- 
tributions of the 7r orbitals, respectively. Panels (a) and (b) 
refer to N2, and panels (c) and (d) correspond to CO, respec- 
tively. The left and right panels display the results using the 
length and velocity form of the dipole operator, respectively. 
The minima due to the two-center interference are indicated 
by green arrows in the figure. 



comparison, see panels (a) and (b) for N2 and panels (c) 
and (d) for CO). 

The contributions from the a orbitals, on the other 
hand, can behave in very different ways depending on 
whether the velocity or the length form is taken, in the 
heteronuclear case, such as CO. For N2, in contrast, dif- 
ferent forms of the dipole operator will mainly lead to a 
different weighting of the contributions from individual 
orbitals to the overall corrections. For instance, in the 
length-form case, the contributions of the 2a u orbital are 
prevalent, while if the velocity form is taken, the contri- 
butions from both a u orbitals are comparable. Another 
interesting feature is that both forms lead to a sharp har- 
monic suppression near 9l = (2n + l)ir/2. This is due to 
the fact that the nodal planes of the \a u and 2a u orbitals 
are parallel to the laser-field polarization for such angles. 

In CO, however, we notice two interesting features. 
Firstly, the nodes disappear if the length form is taken. 
This is partly due to the static dipole moment, which is 
only present in the length form and leads to corrections 
which behave like the HOMO. This implies that the over- 
all dipole will decrease for 6l — it as shown in the figure. 
In contrast, in the velocity form, as this term is absent, 
there are minima near 9l = 7r/4. These minima come 
from the geometry of 2a and 4er, which are the heteronu- 
clear counterparts of \a u and 2er u , leading to a strong 
suppression in the yield close to 8l — (2n+ 1)tt/2. They 
are, however, slightly shifted from these values due to the 



heteronuclear character of the molecule. The remaining 
a orbitals shift these minima even further, to Ql — 7r/4. 



IV. CONCLUSIONS 

In conclusion, our results show that an adequate choice 
of the dipole operator in order to model high-order har- 
monic generation is still an open question. In the single- 
active-electron approximation, enough evidence has been 
provided in [lOJ that the velocity form yields the best 
agreement with the double-slit physical picture. If one 
goes beyond this approximation, however, the velocity 
form leads to a vanishing static dipole moment for het- 
eronuclear molecules. This may be a particularity of the 
basis set involved. 

Another issue in the velocity form are the so-called 
overlap integrals, in which the dipole operator couples 
wavefunctions localized at different ions in the molecule. 
Throughout the literature, it has been argued that such 
integrals are small and may be neglected without loss 
of information. Our results, however, indicate that, de- 
pending on the basis set used, their outcome may be 
comparable to that of the direct integrals, in which wave- 
functions at the same ion are coupled. As the internuclear 
distance increases, however, the overlap integrals become 
vanishing small. 

In contrast, the length- form of the dipole operator 
leads to a non-vanishing static dipole moment for a het- 
eronuclear molecule, and very small overlap integrals, in 
comparison to their direct counterparts. The main prob- 
lem, however, is that such a formulation leads to terms 
in the dipole-matrix elements with increase indefinitely 
with the internuclear separation. They also couple states 
with the same parity, which is unphysical for an odd op- 
erator such as the dipole. Such terms have been also 
found in the single-active electron context [10| , Hence, it 
seems that the velocity form still gives better results. 

Furthermore, suppression of the harmonic spectra 
caused by the presence of nodes in the orbitals is far 
more robust than the maxima and minima coming from 
two-center interference. Interestingly, for heteronuclear 
molecules, despite interference minima being washed out, 
the high-harmonic suppression caused by the geometry of 
the orbitals is still present in the spectra. However, com- 
pared to their homonuclear counterpart, this suppression 
occurs at different alignment angles. 

Finally, in order to incorporate multielectron effects in 
SFA-like models in an appropriate way, one must consider 
the dynamics of the target. Concrete examples are exci- 
tation, relaxation, electron orbits starting or finishing at 
different molecular orbitals, or the motion of the bound 
electrons. This paper has shown that, if multielectron ef- 
fects are considered statically by employing many-body 
perturbation theory around the single-electron strong- 
field approximation, they lead to corrections which are 
several orders of magnitude smaller than the zeroth-order 
spectrum. This holds for both homonuclear and het- 
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eronuclear molecules. In fact, recent results in which ex- where 
citation or relaxation have been incorporated numerically 
exhibited larger corrections [2, @| • 
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APPENDIX 1: OVERLAP INTEGRALS FOR 
DIPOLE MATRIX ELEMENTS 



with 



In this appendix, we are providing the generalized ex- 
pressions for the overlap integrals (|31[) considering the 
dipole operator in its velocity and length forms. 
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A. Velocity form 

We will commence by the velocity-form expressions. If 
only a orbitals are coupled, one may show that 
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If, on the other hand, we consider a transition from a 
a to & tt x or 7Tj, orbital, 
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where C(l a ,l^,l^) = J : (l a )J'(lp,l'y)- The function J-(l a ) is defined in Eq. (IM1) . If only ir x or 7r a orbitals are coupled, 
the overlap integrals read 
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with x = x or y. Finally, if we consider the overlap inte- with 
gral coupling different 7r orbitals, we will find 
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0(J«,J/O = «(Ja,J/»,0,l) (52) 

+H(Z Q + 1, //3, 0, 0) + «(J a , //? + 1, 0, 0) 



and H(l a ,lp, l-y, h) = A(l a ,lp)J(lry, l s ). In Eq. ([52]), A ^ 
x{-2Q'u^G«^/3) + ^(Za,^ - 1,0,0)}. x- The function ^(/ a ,^) is defined in Eq. @. 
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B. Length form 

For the length form of the dipole operator, if only tran- 
sitions involving a orbitals are considered, we find that 

l^'f{a,a) =2^ Jn Jt>\ e 3 '" H l »,h)> (53) 
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For transitions from a er to a 7^ or 7r. y orbital, the overlap 
integral reads 
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with x — X or y. 

Finally, if only n orbitals are involved, one may dis- 
tinguish two types of overlap integrals. Either the dipole 
matrix element couples the same initial and final orbital, 
■k x or 7r y ,or the orbitals tt x and ir\, \ ¥" ^- The former 
case is given by 
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while the latter reads 
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x {£(l a ,lfi)J(0, 0) + ^(J a , lfi)K(0, 0)} , 

with £{l a ,lp) = .4(Z a , ^ + 1) + .4(/ Q + 1, Z^). 

C. Particular integrals J and /C 

Specifically for the basis set employed in this work, we 
find 



AC(1,0) = 



and 
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(66) 

for the integrals K.(l a ,lp). The above-stated expressions 
show that the integrals J~(l a ,lp ± 1) and K.(l a Ap) will 

lead to a polynomial dependence on R/(Q „ +Q> i) 1 ^ 2+ ™- 
This dependence will be damped by the gaussian factor 
in Eq. (|4l)]) , so that the overlap integrals are expected to 
decrease with the internuclear distance. 
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VI. APPENDIX 2: ANGLE-INTEGRATED 
DIPOLE MATRIX ELEMENTS 

In this appendix, we provide the expressions obtained 
for the dipole matrix elements, if the azimuthal angular 
variable is integrated over. Let us start with N2. In this 
case, the HOMO is a 3a g orbital. The orbitals that will 
contribute to the corrections are l7r Ua , , ln u , 2a u and 
\a u . The remaining, a g orbitals do not contribute to the 
corrections for symmetry reasons. We will be able then 
to write for the first-order dipole matrix element, 



dinate </>&, Eq. (|67j) can be written as 

d (1) (k + A(r)) = 2> ff (fc,T,0fc)+2>*.(fc,T,0*)ca#ifc] 
+!>*„(*;, T,0*)sin[0 fc ], (70) 

where V(k 7 T,9 k ) indicates the part of the prefac- 
tors which do not depend on the azimuthal coor- 
dinate cf) k . The modulus square of the product 
[d (1) (k + A(i))] * S 1 ^ (k+ A(i'), when integrated over the 
azimuthal coordinate, leads to the effective prefactor 



d«(k+A(r)) = 4° ff ) o (k+A(r))+dW+4 1 ) +dW , (67) 



where 



v(k,t,t',e k )) = 2nAi(k,t,t',e k ) + 'KMk,t,t',e k )(7i) 
+^A 3 (fc, t, f, e k ) + ^A 4 (fc, t, t', e k ), 



and 



d« = -"£ d 3ag . na ^ nau (k + A(r)) (68) 



d« = -d 3<TgMux ^ ux Qt + A(r)). (69) 



In the above-stated equations, r = t, if and \ — x,y. 
With respect to the dependence on the azimuthal coor- 



where the constants arise from integrals of the form 
J sm m [(/) k ] cos n [4>k]d(j)k, where m,n are even integer 
numbers. The explicit expressions for the functions 
A n (k,t,t',6 k ) read 



Ai(fc,i,t' ',6k) = paa(k,t,t' ,6 k 



(72) 



J 



A 2 (k,t,t',e k ) = \D a ^(k,t,t',e k )\ 2 + \-D^ y (k,t,t',d k )\ 2 + D tTtT (k,t,t',o k )('D* (k,t,t',9 k ) + v* (k,t,t',e k )) 



+v* (T {k,t,t',e k )(v^^(k,t,t',e k )+v 7rv „ v (k,t,t',e k )), 

A 3 (fc, t, t', e k ) = \v^ (k, t, t', e k )\ 2 + \v^ y {k, t, t', 9 k ) I 



(73) 
(74) 



A 4 (fc, t, t', 9 k ) = \V^ y (k, t, t', 9 k ) I + 2?; xTx (k, t, t', 9 k )V ny7Vy (k, t, t', 9 k ) + T> nx7!x (k, t, t', 9 k )V^ yWy (k, t, t', 9 k ), 
and where 



V i:j (k,t,t',9 k ) 



V*(k,t,9 k )V J (k,t',9 k ),i=j 

VX{k,t,e k )Vj{k,if,e k ) + V}[k,t,e k )Vi[k,-e,e k ),i^j. 

I 



(75) 



For CO, there will be two main differences. Firstly, a 
static dipole moment is expected. Secondly, due to the 
fact that the orbitals do not exhibit a well-defined parity, 
more states are coupled by the dipole operator through 
the matrix elements cfoi- In fact, apart from the HOMO, 
there will now be four contributing occupied a states and 
the two degenerate 7r states corresponding to the HOMO- 
1. The explicit expressions for dd will change to 



d£) = - J2 <W.Vw(k + A(r)) + d^, (76) 



with 



d {1) -Vd 

"static 2_^ u na,ii- 



V> 5CT (k + A(r)). 



(77) 



This will lead to changes in T> a (k, r,9 k ). Furthermore, in 
Eq. (f6"9|) 3a g should now be replaced by 5<r and lir ux by 
l7r x . There will be, however, no change in the structure 
of the subsequent equations or in the weights obtained 
upon angular integration. 
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